# Gilardi, Fabrizio, "The Temporary Importance of Role Models for Women's Political Representation", American Journal of Political Science
# Code to replicate Figure 1 (Spatial autocorrelation of the percent of female candidates, 1970-2010)
# gilardi@ipz.uzh.ch, 2014-06-23

# Set working directory
setwd("../Data/")

# Function to convert municipality IDs
bfsnr.convert <- function(i){
	d <- read.csv("nr-bfsnr.csv")
	x <- d$nr[d$bfsnr %in% i]
	return(x)
}

# Load packages
library(spdep)

# Load data
d <- read.csv("dataset-elections.csv")

# Create vector with election years 
yy <- seq(1970, 2010, 4)

# Create matrix to be filled with Moran's I values
moran.cand <- matrix(NA, nrow=length(yy), ncol=4)
rownames(moran.cand) <- yy



## Compute Moran's I for each election year and store it in the matrix created in the previous step

for(s in 1: length(yy)){
	
	# Retrieve election year
	y <- yy[s]
	
	# Convert municipality ID
	nr <- bfsnr.convert(unique(d$bfsnr[d$year==y & is.na(d$pcent.women.cand)==F]))
	
	# Load driving-distance connectivity matrix
	W.auto <- as.matrix(read.csv("W-auto.csv"))
	W.auto <- W.auto[nr,nr]
	
	# Compute inverse such that larger numbers mean closer connection
	W.auto.d <- W.auto
	for(i in 1:length(nr)){
		t <- max(W.auto[i,], na.rm=T)		
		for(j in 1:length(nr)){
			if(W.auto.d[i,j] <= t & W.auto.d[i,j]!=0){ W.auto.d[i,j] <- 1/W.auto.d[i,j] }
			else if(W.auto.d[i,j] > t | W.auto.d[i,j]==0){ W.auto.d[i,j] <- 0}
		}
	}

	# Create vector of percent of female candidates in a given year
	women <- d$pcent.women.cand[d$year==y  & is.na(d$pcent.women.cand)==F]
	
	# Prepare connectivity matrix for Moran's test
	w.n <- mat2listw(W.auto.d)
	listw.n <- nb2listw(neighbours=w.n$neighbours, glist=w.n$weights, style="W")
	
	# Compute Moran's test, extract estimates and store them in matrix
	moran <- moran.test(women, listw=listw.n, randomisation=T, alternative="two.sided")
	moran.est <- moran$estimate[1]
	moran.cand[s,1] <- moran$estimate[1] - moran$estimate[2]
	moran.cand[s,2] <- moran.cand[s,1] - 1.64*sqrt(moran$estimate[3])
	moran.cand[s,3] <- moran.cand[s,1] + 1.64*sqrt(moran$estimate[3])
	moran.cand[s,4] <- moran$p.value
	
}


# Plot Moran's I

#pdf(file="Figure-1.pdf", paper="special", width=7.5, height=5.5)
par(mar=c(3.5,3.5,2,1), mgp=c(2.5,0.8,0), cex.axis=1.2, cex.lab=1.2, font.main=1)
plot(c(1970, 2010), range(moran.cand[,2:3]), type="n", axes=F, ylab="Spatial autocorrelation (Moran's I)", xlab="")
axis(1, at=yy)
axis(2)
points(yy, moran.cand[,1], type="b", pch=1, cex=1.5)
segments(yy, moran.cand[,2], yy, moran.cand[,3])
abline(h=0, lty=2)
#dev.off()
